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ABSTRACT 

A finite difference procedure for the computation of stress and deformation in a 
three-dimensional elastic-plastic solid is presented. This iterative technique, based on 
Newton’s method for determining the roots of a system of polynomials, is applicable to 
both linear and nonlinear problems. The procedure requires a minimum of computer 
storage capability and moderate computer running time. 

Three sample stress -concentration problems are treated. The first problem is the 
analysis of a plate in plane strain with a slit at its center, the second is the analysis of a 
thick plate with a rectangular slit through its thickness, and the last is the analysis of a 
thick plate containing a semielliptical slit which extends only halfway through the thick- 
ness. The loading condition in all sample problems is uniaxial tension normal to the 
slit. 
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A NUMERICAL PROCEDURE FOR CALCULATING STRESS AND DEFORMATION 
NEAR A SLIT IN A THREE-DIMENSIONAL ELASTIC- PLASTIC SOLID 

by David J. Ayres 
Lewis Research Center 

SUMMARY 

A finite difference procedure for the computation of stress and deformation in a 
three-dimensional elastic-plastic solid is presented. This iterative technique, based on 
Newton’s method for determining the roots of a system of polynomials, is applicable to 
both linear and nonlinear problems. The procedure reqmres a minimum of computer 
storage capability and moderate computer running time. 

Three sample stress -concentration problems are treated. The first problem is the 
analysis of a plate in plane strain with a slit at its center, the second is the analysis of a 
thick plate with a rectangular slit through its thickness, and the last is the analysis of a 
thick plate containing a semielliptical slit which extends only halfway through the thick- 
ness. The loading condition in all sample problems is uniaxial tension normal to the 
slit. 


INTRODUCTION 

Numerical methods which involve the solution of large systems of simultaneous 
equations have long been applied to two-dimensional problems in solid mechanics. The 
relaxation technique of Southwell (ref . 1) , the matrix displacement methods of Clough 
(ref. 2) and Argyris (ref . 3), the discrete model of Harper and Ang (ref. 4), and the 
method of successive elastic solutions of Mendelson (ref. 5) are some of the procedures 
that have been successfully applied to two-dimensional elastic -plastic problems. 

The matrix displacement method is the only numerical procedure that has been suc- 
cessfully applied to three-dimensional problems (ref . 3). The computer storage and 
running time requirements, however, are great because of the large number of simul- 
taneous equations to be solved. In particular, the accurate determination of the stress 
near a slit in a three-dimensional elastic-plastic body requires the solution of so many 



simultaneous equations that the running time on present day computers discourages the 
use of this method. 

In order to efficiently solve three-dimensional problems using computers that are 
available today, a new procedure is required. This procedure must permit the solution 
of thousands of simultaneous equations fairly rapidly. It must be sufficiently general so 
that problems concerning many different types of materials can be solved, and it must be 
sufficiently accurate to give a useful approximation of the stress close to a slit. This 
report presents such a numerical procedure. 

The procedure presented herein is an extension of the general method for computing 
stress and deformation developed in reference 6 and was inspired by the iterative method 
of Harper and Ang (ref. 4). For this procedure the partial differential eqiiations of equi- 
librium of the solid are written as finite difference equations in terms of the displacement 
of points of the solid. These equations are solved iteratively by a relaxation technique. 
This relaxation converges for linear and nonlinear equations in three dimensions if a 
reasonable first approximation to the displacement of the solid is assumed. The stress 
in the body is computed from the resulting displacements. 

In this report, the solid is assumed to be elastic and perfectly plastic, obeying 
Hookes’ Law, the Mises-Hencky yield condition, and the Prandtl-Reuss flow rule. The 
procedure is illustrated by three examples of a slit in an elastic-plastic plate. The first 
example is a slit in plane strain, the second is a thick plate with a rectangular slit 
through the thickness, and the third is a thick plate with a semielliptical slit only par- 
tially through its thickness . 


SYMBOLS 

In this report the conventional summation notation is employed. For example, 

^ii " ^11 * ^22 ^ ^33 

where the range of the index is three unless otherwise indicated. The following symbols 
are used: 

b width of plate 

c half length of slit in the Xj direction 

d small increase in a displacement component 

E Young's modulus 

EQ equilibrium equation 


2 



e 


volumetric strain 


G shear modulus, G = - — — — 

2(1 + v) 

i, j indices 

Jg second deviator stress invarient 

E 

K bulk modulus, K = 

3(1 - 2v) 

yield condition, k^ = i cr^ 

3 ^ 

m, n index limits 

p grid point 

R residual (amount by which a finite difference equation is not satisfied) 

RF relaxation factor 

S hydrostatic stress 

T stress tensor 

TB traction boundary equation 

U displacement 

W plastic work 

X distance 

A difference between values at two sucessive load levels 

6.. Kronecker delta, 6.. = 

e strain 

V Poisson's ratio 

CTy uniaxial yield stress 


1 if i = j 
0 if 1/j 
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BASIC EQUATIONS 


Every point in a solid body must satisfy the equations of equilibrium. For a station- 
ary body with no body forces acting these equations can be written in Cartesian coordi- 
nates as follows: 


3T. 


ax. 


ii = o 


r= 1,2,3 
1= 1,2,3 

V. 


( 1 ) 


where the T. . are the components of the stress tensor T in the body at the point with 
coordinates (Xj^, X2, Xg) . 

The behavior of each material can be described by a constitutive equation. An elas- 
tic, perfectly plastic material is described by the combination of Hooke’s Law, the 
Mises-Hencky yield condition, and the Prandtl-Reuss flow rule. In Cartesian coordi- 
nates the constitutive equation derived in reference 4 can be written in incremental form 
as 


where 


ATy = 2G 


Ae.. - 6.. Ae - 
1] 1] 


AW 


2k 


(T.. - 6..S) 

2 13 13 


+ 3K5.. Ae 


e = 


^ii 


s = 


Tii 


AW = 


(T.. - 6..S)(Ac.. - 6.. Ae) if Jg = k^ 


0 if Jg < k^ 


( 2 ) 


^ 2 %(Ti 3 


6..S)(T.. - 5..S) 
1 ] 1 ] 13 ' 


( 3 ) 




The Tj. are known from the previous load increment, and U- are the components of the 
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displacement of the body at the point (X^^, Xg, Xg). On the boundary of the body, either 
traction or displacement must be assigned. 

If the stress tensor T satisfies equation ( 1) for a particular boundary condition and 
constitutive equation and the boundary condition is altered so that the stress changes by 
AT, from equation (1) the AT.j must satisfy the equation 




(4) 


Substituting equations (2) and (3) into equation (4) results in a set of three simulta- 
neous partial differential equations in terms of the three unknown displacement compo- 
nents and the known stresses and displacements from the previous load increment. 

These equations are solved by a finite difference technique. Certain points in the body 
are designated as grid points, and the displacement U of the body is defined at these 
points. When the finite difference approximations to the derivatives of U are substi- 
tuted into equation (4), written at each grid point in the body, the system of partial dif- 
ferential equations reduces to a large number of simultaneous polynomial equations which 
are written symbolically as 

EQip(AUi, AUg, AUg) = 0 

where n is the number of interior grid points. 

The constitutive equation must be satisfied at every boundary grid point. When 
traction is prescribed in a particular direction at a boundary point, the constitutive equa- 
tion is used to determine the displacement in that direction. When the finite difference 
approximations to the derivatives of U are substituted into equation (2), the traction 
boundary condition is written in symbolic form as 


fi = 1,2,3 
|p = 1, . . . , n 


i = 1,2,3 

TB. (AU., AUg, AUg) = 0 < (6) 

^ ^ |p = 1, . . . , m 

where m is the number of boundary points where traction is prescribed. The displace- 
ment U at each grid point is determined by the simultaneous solutions of equations (5) 
and (6) . 
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NUMERICAL PROCEDURE 


The numerical procedure for the solution of the finite difference equations (eqs. (5) 
and (6)) is briefly stated as follows: 

(1) Designate certain points of the body as grid points. Choose a set of values for 
the displacement U that approximates the effect of the applied load at all these points. 

(2) Starting at a convenient point (preferably near the slit), substitute the displace- 
ment components into the equilibrium equation (eq. (5)) or the boundary condition (eq. (6)) 
for a particular direction i. Then calculate the residual R. ^ , the amount by which the 
equation is not satisfied. 

(3) Alter the displacement component in the i direction U. by a small arbitrary 
amount d. Compute the residual R. ,, when the altered displacement is used. 

(4) Calculate a new value that causes the residual for this equation to be less 

than R. ■. found in step 2. This value is calculated by the extrapolation illustrated in 
figure 1 to be 


u = u + = u 

i(new) i(guess) i(guess) 


R. 






( 7 ) 


(The selection of the relaxation factor RF is discussed for each example. ) 

(5) Repeat this procedure (steps 2 to 4) for each direction i at each point p in the 
body. 

(6) Repeat the complete procedure (steps 2 to 5) until all the residuals are small. 

The procedure is said to be converging when the sum of the absolute values of the resid- 
uals is becoming smaller. It is assumed to have converged when this sum becomes less 
than a prescribed small number. 

This procedure is a form of Newton’s method (ref. 7) for determining the roots of a 
system of polynomials and, therefore, does not reqvdre the constitutive equations to be 
linear in displacement derivatives. In this report, however, only the linear equation (2) 
is used. The one requirement for the procedure is a reasonable first guess of the dis- 
placement of the body due to the boundary traction. It is shown in the examples, however, 
that these assumed values need not be very close to the true values. 
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EXAMPLES 


The numerical procedure is used to obtain approximate solutions for three stress- 
concentration problems. The first example illustrates the two-dimensional plane strain 
analysis of a rectangular plate which contains a central slit, the second illustrates the 
more practical three-dimensional analysis of a plate of finite thickness containing a rec- 
tangular slit through its thickness, and the third illustrates the case of a thick plate con- 
taining a semielliptical slit which extends only halfway through the thickness. The di- 
mensions of these three plates are shown in figure 2. The elastic-plastic material prop- 


erties for all examples are as follows: 

Young’s modulus, E, ksi (N/m^) .................... 3x10^(2.06x10^^) 

Poisson’s ratio, y ................................. . 0. 3 

Uniaxial yield stress, Oy, ksi (N/m^) 200(1,38x10^) 


Plane Strain 

The geometry of the first example is illustrated in figure 2(a), A uniform tension is 
applied in the X£ direction. No deformation occurs normal to the plane of the plate. 

The slit half-length c equals 1/6 of the plate width b. 

Figure 3 shows a composite of the three successive grid patterns used to compute 
the elastic solution on one-quarter of the symmetrical plate. The method of choosing the 
successive grids is briefly stated as follows: 

(1) The entire quarter plate is divided into squares with side length b/18 (grid 1) . 

The uniform tensile load is applied to the upper boundary of the plate (boundary line 1). 

In this example the displacement values for the plate without a slit are assumed for the 
first trial. The elastic solution for this grid and boundary condition is determined within 
a certain error by the numerical procedure. 

(2) Stress values calculated on grid 1 are interpolated along boundary line 2. The 
displacements on grid 1 are interpolated to give an approximation to the solution of the 
elastic plate below boundary line 2 on the square grid with side length b/36 (grid 2). 

The elastic solution of this separate problem with a traction boundary condition pre- 
scribed at boundary line 2 is computed. 

(3) Stresses calculated on grid 2 are interpolated along boundary line 3. The dis- 
placements on grid 2 are interpolated to give an approximation to the solution of the elas- 
tic plate below boundary line 3 on the square grid with side length b/72 (grid 3). The 
solution of this final problem is considered the elastic solution for the plate in plane 
strain. 
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The technique of placing a more refined grid over a smaller area could be continued 
many times. These three grid sizes, however, establish the trend to be expected if 
more refined grids are used. The technique permits the use of a fine grid where the 
stress is expected to vary greatly with distance and a coarse grid wher e the stress is ex- 
pected to be more nearly constant. 

The stress ahead of the slit and normal to it calculated on each of the three square 
grids is presented in figure 4 and table I. The solutions of the elastic plane problems 
for an infinite plate with a slit b/3 in length (ref. 8), an infinite plate with a series of 
colinear slits b/3 in length with centers b apart (ref. 9), and an unpublished solution 
for the finite-width plate by A. Mendelson of the Lewis Research Center are shown for 
comparison. The variation of the stress distribution with grid size sets an obvious trend 
which appears reasonable when the analytical solutions are considered. 

In order to calculate the analytical solutions, the exact location of the end of the slit 
must be known. The numerical procedure requires only that the grid points on the slit 
be considered traction-free boundary points and that the grid points not slit be considered 
interior points. The end of the slit then lies somewhere between two grid points. The 
xmiform tensile load deforms the slit into a nearly elliptical shape. Extrapolation of the 
displacement of the grid points on the slit provides an estimate of the location of the tip 
of the slit where no displacement in the Xg direction occurs. This value is approxi- 
mately one -quarter of a grid space before the first interior grid point in the line of the 
slit. This more precise crack length of (11. 75/12. 00) (b/3) is used to calculate the ana- 
lytical solutions shown in figure 4. 

The first grid point (on grid 3) yields at an applied gross tension of 0. 424 cr^. The 
plate above boundary line 3 is assumed to remain elastic and to remain unaffected by the 
growth of the plastic zone near the tip of the slit. After initial yielding, the load on 
grid 3 was increased by 3 percent of the previous load for 28 increments. In order to 
obtain a fairly rapid solution, a fixed small number of iterations (30) was employed for 
each increment. The growth of the plastic zone is illustrated in figure 5 for various ap- 
plied tensile loads. The redistribution of stress ahead of the slit and normal to it due to 
the plasticity of the material is illustrated in figure 6 for one particular load. 

The convergence of the elastic solution on the three successive grids is shown in 
figure 7. In this figure the mean residual (i. e. , the sum of the absolute value of the re- 
siduals divided by the number of grid points) is plotted against increasing number of it- 
erations. For each grid the mean residual of the finite difference equations decreases 
with increasing iterations. As the grid boundary comes close to the slit, the mean re- 
sidual increases because a greater percentage of grid points are very near the slit. 

The relaxation factor chosen for this two-dimensional problem was 1. 0 because 
numbers greater than 1.0 decreased the rate of convergence and numbers less than 1. 0 
caused oscillation or divergence of the solution. The iteration was stopped for each grid 
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when the increment AUj to every displacement component U. satisfied the condition 


Ui fi 

AU. < - — — + 10“® 

^ 4 

lx 10^ 


( 8 ) 


-fi 

The second term (10" ) on the right side of expression ( 8 ) was chosen to be about the 
same magnitude as the first term. This expedient is needed for those stations where the 
displacement is zero. 

Figure 8 illustrates the growth of the mean residual with the number of load incre- 
ments. As the plastic zone size increases, the mean residual becomes large because a 
fixed number of iterations are employed. The large zones illustrated in figure 5 are cal- 
culated from stress values that may contain five times as much error as the elastic solu- 
tion. The residuals could be reduced to any prescribed amount, however, by allowing a 
greater number of iterations for each loading increment. The largest zone, at a load of 
0. 970 Oy, reaches boundary line 3 and, therefore, violates the assumption that boundary 
line 3 would remain unaffected. These large zone sizes then must be considered only as 
trends. 


Rectangular Slit 

The geometry of the second example is illustrated in figure 2(b). The thick plate 
contains a rectangular slit of length b/3 which extends entirely through its thickness. 
Figure 9 shows, on one-eighth of the symmetrical plate, the three successive grids that 
were employed. The technique of refining the grid near the slit is exactly as discussed 
for the plane strain example. In this case the extent of each grid was determined pri- 
marily by computer core storage limitations. 

The elastic stress ahead of the slit and normal to it due to the uniform tensile load 
in the Xg direction is shown in figure 10 for the center of the plate (Xg = 0) and for the 
face (Xg = b/12) for each of the three successive grids. The stresses show the same 
trends with successively smaller grid sizes as those in the plane strain example. The 
variation of the stresses T^j^, Tgg, T 33 through the thickness of the plate is illustrated 
in figure 11 and compared with the plane strain values calculated in the previous ex- 
ample. These results demonstrate the same trends as the stresses near a hole in a 
thick plate calculated by other authors (refs. 10 and 11 ). 

The first grid point (grid 3) yields at an applied gross tension of 0. 338 Uy. Since the 
boundary line 3 is near the slit, it is not reasonable to assume that the material on that 
line would be unaffected by the growth of the plastic zone. Therefore, it is assmned that 
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the plate above boundary line 2 remains elastic and unaffected by the plastic zone. For 
each load increment then, the elastic-plastic solution on both grids 2 and 3 must be com- 
puted. After initial yielding, the uniform tension was increased by 3 percent of the pre- 
vious load for 30 increments. The growth of the plastic zone normal to the slit at the 
center of the plate (Xg = 0) and at the face (Xg = b/ 12) is shown in figures 12 and 13, re- 
spectively. 

Figure 14 shows the growth of the plastic zone ahead of the slit in the plane of the 
slit. This figure shows that the plastic zone does not extend to the center of the plate 
when the zone is small. At larger applied loads, the zone extends through the plate and 
becomes greater at the centerline than at the faces. 

The convergence of this example is similar to the convergence of the plane strain 
example. Expression (8) is employed as the convergence criterion for the elastic solu- 
tion. Thirty iterations are performed for each load increment; therefore, the growth of 
the mean residual with increasing plastic zone size is also similar to that described in 
the previous example. 

A relaxation factor of 1. 4 was chosen to assure convergence of the procedure. For 
three-dimensional problems the relaxation factor 1. 0 causes rapid divergence. Relax- 
ation factors larger than 1. 4 assure convergence at a slower rate. 


Semielliptical Slit 

The geometry of the final example is illustrated in figure 2(c). The thick plate con- 
tains a semielliptical slit with a major half axis of b/6 and a minor half axis equal to 
b/12 or half the thickness of the plate. As in the other examples, three successive grids 
were employed to compute the elastic stress distribution. The extent of the grids, de- 
termined primarily by computer core storage limitations, is shown on one-quarter of the 
symmetrical plate, in figure 15. The initial trial displacements on grid 1 are approxi- 
mations to the analogous plane stress problem in each (X^, Xg) grid plane in the and 
Xg directions and zero in the Xg direction. 

Contours of elastic stress in the plane of the slit and normal to it due to the uniform 
tensile load in the Xg direction are drawn in figure 16. The first grid point (grid 3) 
yields at an applied gross tension of 0. 437 Oy. As in the previous example, the part of 
the plate above boundary line 2 is assumed to be unaffected by the plastic zone growth, 
and the elastic -plastic solution must be computed on both grid 2 and grid 3 for each load 
increment. After initial yielding, the uniform tensile load is increased by 3 percent of 
the previous load for 24 increments. The plastic zone growth on the face containing the 
slit is illustrated in figure 17. The zone growth through the thickness to the back face of 
the plate is shown in figure 18. Notice that the zone grows from the slit tip and then 
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joins with another zone initiated on the back face. The stress on the back face directly 
ahead of the slit is still at a low elastic stress. A three-dimensional sketch of the plastic 
zone at applied stress of 0. 888 ov is shown in figure 19 for one-quarter of the symmet- 
rical plate. 

The convergence of this example required a relaxation factor of at least 1. 4. As in 
the previous three-dimensional example, values less than 1. 4 caused rapid divergence. 
Expression (8) was employed as the convergence criterion for the elastic solution, and 
30 iterations were performed for every load increment beyond the initial yield point. 

The mean residual thus behaves similarly to that discussed in the plane strain example. 
The mean residual growth during increasing load and plastic zone size is illustrated in 
figure 20. As discussed in the section Plane Strain, a greater number of iterations for 
each increment could reduce the mean residual to any desired amount. 


DISCUSSION 

There are several requirements that a new numerical procedure must fulfill. First, 
it must produce an answer, that is, it must converge. In the three examples illustrated, 
the convergence is monotonic; thus, the error is controlled by the number of iterations. 

A relaxation factor, based on the number of dimensions, grid element shape, Poisson’s 
ratio, and perhaps other parameters, must be chosen wisely. In all cases attempted, 
relaxation factors of 1. 0 and 1. 4 appear satisfactory for two- and three-dimensional ex- 
amples, respectively. These values were determined experimentally, however, and no 
theoretical justification has been made. The requirement that a reasonable first approxi- 
mation be made is apparently not harsh since the trivial value of uniform tension permits 
convergence. 

Another requirement of a numerical procedure is that it produce results of useful 
accuracy. The present numerical results can be compared to analytical solutions for 
plane examples, and the trends in some three-dimensional problems can be compared 
with other problems that have been solved. The numerical results of the plane strain 
example compare favorably with the analytical solutions of similar problems. The stress 
variation through the thickness of the plate with the rectangular slit illustrates the same 
trends that occur in the exact solution of other related problems. 

A final criterion by which to judge a numerical method is economy. Table II lists 
the computer running time on the NASA Lewis IBM 7044/709411 for the three examples. 
These times are much less than those required to solve problems involving fewer un- 
knowns by the matrix displacement method of Argyris (ref. 3) . In order to obtain better 
precision, however, times much greater than those listed may be required. Since it is 
less expensive to use a smaller computer, storage minimization is also an economy. 
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The only quantities that need to be stored for each grid point at any particular load incre- 
ment are six stress and three displacement components calculated for the previous in- 
crement and three displacement and three flexibility coefficient components (the quanti- 
ties (R. , - R. o)/d) for the current increment. The most important economy is conven- 
ience, however, since the method can be employed on any farily large computer with no 
peripheral storage or other special equipment. Although the problems considered herein 
have not included strain hardening, the method is sufficiently general to allow the use of 
various material properties including arbitrary strain hardening. 


CONCLUDING REMARKS 

The iterative procedure presented herein was found to be convergent. The error in 
the solution of the difference equations is known after every iteration and can be con- 
trolled by the number of iterations. The procedure can be i:^ed on most large computers 
presently available because a minimum of core storage and no peripheral storage is re- 
quired. The running time for problems involving thousands of simultaneous equations is 
moderate. The stress values calculated in the elastic plane strain analysis of a plate 
containing a central slit were compared with an analytical solution of Westergaard and 
are in good agreement. Although no exact comparisons are available for the three- 
dimensional examples, high confidence is placed in the usefulness of the procedure be- 
cause of its performance on the two-dimensional problem. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, April 23, 1968, 

124-08-06-01-22. 
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Figure 3, - Grid on one-quarter of plane strain example. 
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Figure 4. - Ratio of elastic stress ahead of slit in direction of load to remotely applied stress for plane strain 
example. 
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Figure 5. - Plastic zone size for plane strain example. 
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Figure 6. ~ Ratio of elastic-plastic stress ahead of slit in 
direction of load to remotely applied stress for plane strain 
example. 
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Figure 16. - Elastic stress in direction of load in plane of semielliptica! slit. Contours of ratio of stress 
to remotely applied stress. 
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Figure 20. - Ratio of accumulated mean residual to 
mean residual for elastic solution for grid 2 of elastic- 
plastic semieMiptical slit example. 
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